Majorana path integral for nonequilibrium dynamics of two-level systems 
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We present a new field-theoretic approach to anaylize non-equilirbium dynamics of two-level 
systems (TLS), which is based on a correspondence between a driven TLS and a Majorana fermion 
field theory coupled to bosonic fields. This approach allows us to calculate analytically properties 
of non-linear TLS dynamics with an arbitrary accuracy. We apply our method to analyze specific 
TLS dynamics under a monochromatic periodic drive that is relevant to the problem of decoherence 
in Josephson junction qubits. It is demonstrated that the method gives the precise positions of the 
resonance peaks in the non-linear dielectric response function that are in agreement with numerical 
simulations. 
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I. INTRODUCTION 

A driven two-level-system (TLS) represents a canoni- 
cal dynamical system that features a rich variety of in- 
teresting non-linear phenomena including various reso- 
nance effects, quantum interference phenomena due to 
level crossings, coherent destruction of tunnelling, etc 
(see e.g., Ref. [l|). Despite the formal simplicity of its 
formulation, this quantum dynamical system does not 
admit an exact solution in a closed analytical form for 
an arbitrary external drive and one often resorts to nu- 
merical simulations or various approximation schemes 
for its analysis. There are deep mathematical reasons 
for the lack of our ability to solve the corresponding 
differential equations, which go back to the old works 
on Riccatti differential equations and Lie theory. The 
mathematical problem itself has a wide spectrum of ap- 
plications in technology and physics ranging from tech- 
nologically important nuclear magnetic resonance spec- 
troscopy^ to Maxwell-Bloch theory of two-level lasers^, 
non-equilibrium superconductivity^, and fundamental 
field-theoretical models such as the Wess-Zumino-Witten 
theoryS.. Furthermore, with the increased technological 
ability to fabricate and control quantum systems, new 
realizations of the model arise on a continuous basis, 
such as, for example, artificial "atoms" interacting with 
strongly oscillating fields^.. This variety of new applica- 
tions has motivated focused theoretical researches of the 
model recently, see e.g., Refs. [sL^d]! that studied a weak 



driving limit with small number of photons and Ref. 12 1 
that investigated the regime of strong driving with large 
photon numbers, just to name a few relevant papers. 

Another important field, where the problem of driven 
TLS dynamics is of great importance, is quantum com- 
puting. There are various realizations of qubitai^^— , 
which in the course of quantum evolution may exhibit 
interesting non-linear phenomena, such as interference 
between multiple Landau-Zener transitions at a level 
crossing, where adiabatic evolution between them results 
in an oscillatory qubit magnetization in the regime of 
strong qubit drivin g^^d"^ . Physics of driven TLS shows up 



in qubits also from a different perspective: Low-energy 
charge defects are widely believed to be the dominant 
source of dephasing in superconducting Josephson junc- 
tion qubitai^. There, interactions between the charged 
TLS defect and an applied electric field gives rise to the 
same problem of a non-equilibrium TLS under a periodic 
monochromatic drive. In general driven TLS was a sub- 
ject of intensive studies during the last decade, different 
aspects of which are presented in e.g., Refs. [l9l - l23l 

In most cases mentioned above, the basic problem that 
we are actually interested in is summarized by the simple 
Hamiltonian, H = — A(t) -o-, which leads to the following 
evolution operator U{t) and the "partition function" Z 



U{t) = Texp 



dtA{t) 



Z = tr[/(T). (1) 



Here T is the time ordering operator, & — (cti, (T2, o's) 
is the vector of Pauli matrices, and A(t) — 
[Ai(t), A2(t), A3( t)] is the three-component external 
driving field. 

In this paper, we propose an analytical Majorana field- 
theory approach, which is capable of describing arbitrary 
external driving fields. Though the approach can be 
applied to the most general case however we consider 
following form of the periodic driving field: A(t) = 
[Ai(t), A2(t), A3], where A3 is time independent con- 
stant (this assumption does not break the generality as 
it corresponds to transforming to a frame with static 
A3) and we also assume that the two time- varying com- 
ponents have a finite number of harmonics, Ai{t) = 
Ci + cos{ujnt), i = 1,2. Lct us first summarize 

three steps that one has to follow to analyze this quan- 
tum dynamical system within our approach: (1) Fourier 
transform the fields: 



Ai{uj) = ei(5(iLj) + 



E 

n,fT— lb 



An6{uj + aujn); 



(2) 



(2) Write exact Dyson equations for the fermionic 
Green's function, /C(a;,a;') in a Majorana field theory; 
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(3) Solve the resulting equations recursively. The follow- 
ing presentation deciphers this prescription and provides 
a specific example of its use for the most experimentally 
relevant (but analytically unsolvable) case of a simple 
monochromatic drive: 



e A 
^iW = ^+Acos{wt), A2{t) EE 0, Asit) = — = const, 



(3) 



were e is an energy splitting, w and A are the frequency 
and amplitude of an external field. Our main result is 
a practically-useful continued fraction representation of 
the correlation function, K{t,t') — tr((T3(t)(T3(i')), that 
contains its exact spectrum. 



Equivalence of the driven TLS ([3]) and the Majorana field 
theory ^ can be shown in three steps: 

(i) First we expand the exponent of the interaction 
term in the expression of the partition function, Z [see 
Eq. Q] 



dt„ 



{Ap{tn)^p{tn)Utn) ~ ^C^.^M (^n ) A3) 



X exp 



dt 



(6) 



II. ACTION IN TERMS OF MAJORANA 
FERMIONS 



(ii) At the second stage we make use of the free Green's 
function of the Majorana fields (see details in Appendix 
A) 



Our work is based on the observation that the ex- 
pression for the "partition function," Z, resembles the 
Green's function of a spinning particle passing in a one- 
dimensional space from point x{t = 0) to x(t = r) in 
Feynnian path integral quantization approach j24| , where 
A(t) has a meaning of the velocity, x(i), of the spin- 
ning particle. Subsequently, the conversion to Majorana 
fields is achieved by using the approach of Refs. 24,25], 
which provides a simple prescription: the Pauli matrices 
o'fj,, ^3 — should be replaced by the Majorana fields 

^fj,{t), ^3{t) respectively. In other words, Pauli matrices 
can be regarded as a quantized version of path-integral 
Majorana fields. We note that Majorana fermion is its 
own antiparticle, i.e., its creation and annihilation oper- 
ators are identical. As fields, they can be described by 
real valued Grassmann variables^^*^ denoted below as 
^(t). The spin dynamics has been investigated in a dif- 
ferent context using Majorana fermion representation in 
Ref. '2?]. Similar approaches have been employed previ- 
ously in Refs. |28l - [34j . 

A cornerstone of this work is that the "partition func- 
tion," Z, i.e., the trace of the evolution operator can 
be exactly reproduced within a field theory of three Ma- 
jorana fermions ^i(i), ^(i), and £,3{t), defined by the 
functional integral 



(4) 



where Zq = Det [jj] , and the action [here we set A3(t) 
A/2 = const] has the form: 



«5({a,A) = 



dt 



1 



+ A^(oe^(06(0 - -tMt)Ut)^ 



.(5) 



In Eq. © fi,u = 1,2, ei2 = -621 = 1, en = €22 = 
and the summation over repeated indices is implied. 



(7) 



X exp 



= Sp^sign[ti - t2], 



This form of the Green's function follows straightfor- 
wardly from the identity -^sign[t — t'] ~ d[t ~ t']. Then, 
by use of the Wick's theorem, one can replace correla- 
tion functions of the Majorana fields in the expression 
([6|) by the trace of the product of Pauli matrices. This 
procedure yields 



X exp 



dt 



allpairings 

= Tr[ap^ap^ ■■■crp„], 



(8) 



where P is the number of permutations needed for ob- 
taining the set of indices //i/iii/i2Mi2 ' • ' MTi/2/^i„/2 from 
/ii/i2 ■ ■ ■ fJ-n (note that n here is even). It is also impor- 
tant to note that sign[t — t'] factors in the Majorana 
fermion Green's function ([7]) disappear in the expression 
above due to the time ordering of fields in Eq. © . Simi- 
lar situation is with the functional integral over the fields 
^3, which satisfy the condition (6(06(^')) = sign[t~t']. 
This suggests that the fields £,3{t) simply can be effec- 
tively dropped from Eq. (|5]), as the whole integral gives 
one. Finally, the expression ([8]) shows that one can re- 
place Majorana fields, ^p{t), by the corresponding Pauli 
matrices, ct^. 

(iii) Upon replacing the fields £,p{t) in Eq. ^ by the 
Pauli matrices, cr^, dropping ^3(1) and subsequently col- 
lecting the obtained series back to the exponent one will 
obtain Z = TrTexp {i dtA{t) • o-}. In derivation of 
this formula we also use the relation (T3 = iaia2. 
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As we mentioned above, action ([5]) can be interpreted 
as the path-integral quantization of relativistic spinning 
particle in a fixed gauge in an external magnetic field 
A3(i) = A/2^. 

Finally, Gaussian integration over the fields 
/i = 1,2 gives the partition function in the form 2, — 
J D[£,3it)]exp{iS3[(3{t),A{t)]}, where the effective ac- 
tion for the £_3{t) field reads 



nian reduces to \3l 



iSs Ks, A(t)] = - / dt^3C3 + 2 / dtdt'Ut') X (9) 
4 Jo Jo 



Here G±{t - t') = \e^^^<^-^'hig-a{t - t'] is the Green's 
function of the differential operator (9t±iA), while 
A±(i) = Ai(t)±?A2(t). Eq. © is one of the new results 
presented here. 

For our purposes it is convenient to express 1S3 through 
the Fourier images of the fields and functions in Eq. ([9]) . 
Then, in the limit t — > cx), we obtain for 53 



53K3,A(t)] = -^ / du:dLo'^3H1C~\uj,uj')^3{^'), 



(10) 



where IC{uj,uj') — [ll!6{uj + uj') + G{ll!,ll!')] ^, and 
G{uj,uj') is an antisymmetric kernel given by 



Zbcs = / 2?[A,A,c,c+]e-^-^^ 



Sbcs ^ / dr|^V^k(5. + /ik)^k + -AA| (12) 



where g is the interaction constant, 

■0k = 

defines the Nambu spinor and 

/ik = CkCTS + AiCTi -I- A20-2 




(13) 



(14) 



is the matrix Hamiltonian, where A = Ai — iA2, A = 
Ai -f iA2. Remarkably, the form of the operator /ik is 
very much reminiscent to our driven Hamiltonian ([3]) , but 
with the third component, ek, being a time independent 
constant. Therefore it is straightforward to employ the 
above developed technique of Majorana field theory to 
represent the exact BCS partition function, Zbcs as a 
functional with respect to one specie Majorana fermion, 
^3(i), and Hubbard- Stratonovich boson fields, A and A, 
see Eq. (fTO)) . If one is interested for example in calcula- 
tion of Zbcs J then, by integrating out £,3{t) one produces 
an effective bosonic action: 



G{uj,uj') = 2 doj 



A+(wi)A_(w +LU1+ lu'){lu' - LU) 



(uji + UJ + A){uji + A) 



(11) 



Due to the novelty of our approach, we first briefly de- 
rive several established results and then turn to the main 
problem of dephasing in the presence of a monochromatic 
drive. 



III. APPLICATION TO 
SUPERCONDUCTIVITY AND VERIFICATION 
OF THE BCS RESULT 



The aim of the present section is to apply the devel- 
oped technique to study the nonperturbative properties 
of the pairing Hamiltonian. The textbook expression of 
the partition function of the pairing model is given in 
terms of a functional integral with respect to complex 
Grasmann fields, Ckcr(i), Ckcr(i), cr =tj4'i with the action 
including four-fermionic pairing interaction. The stan- 
dard approach to treat this action is to decouple the 
interaction term by introducing a set of bosonic fields, 

A{t), A{t), over which one will have an additional func- 
tional integral. Partition function corresponding to the 
zero-particle and paired sectors of the pairing Hamilto- 



Zbcs^ J VA^e^^ff<^^^ 
_A2(e) 1 



(15) 



+ - tr log [ojS{uj + Lj') + G{uj, uj')] 
9 2 



where G{uj,uj') is a functional of A{t) and A{t), and is 
defined by Eq. ([TT|) . Note that the functional integral 
Eq. (jlSp is formally nothing but a nonlinear functional 
determinant written in energy (rather than imaginary 
time) space. Variation of the effective action with re- 
spect to the bosonic fields gives the equation of motion, 
^^^' f — 0, which in turn yields the general gap equation 



+ i / dujduj'IC{uj,uj') 
9 2 7 

A±{Cj + oj + uj'){oj' - oj) 



(16) 



[UJ - 



ek)('^ 



^k) 



written in the energy representation. 

The pairing hamiltonian itself is designed to describe 
the superconducting phase, where the order parameter, 
A, is different from zero. The gap equation (|16p . i.e. the 
equation for the Fourier-transformed bosonic Hubbard- 
Stratonovich field, A{uj), describes the dynamics of the 
order parameter in the global gauge symmetry broken 
phase. 
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The BCS mean- field solution corresponds to the choice 
A±(u!) — Aq (5(w), where is uj independent. In other 
words it assumes that the time-dependence of the order 
parameter is unimportant. Then, from Eq. (Ilip we find 
that 



O 



""^o Siu; + oj'). (17) 



and 



A„MFiW,W ) — 2 5 



S{uj + uj'). (18) 



Making use of Eqs. (|17p and (fT8)) . and substituting them 
into (ITB|) . we obtain 

--+iV(0) /do.^— 1 ^ = 0, (19) 

which reproduces the standard BCS gap equation with 
A^(0) being the approximated to a constant density of 
states. 



IV. NON-DISSIPATIVE TWO-LEVEL SYSTEMS 

Neglecting environmental effects on the driven TLS we 
consider the Hamiltonian H = — A(t)o-, with the driving 
fields given by Eq. ([3]). Our goal is the calculation of 
spin-spin correlation function 

K{t,t') = tT[U{-oo,t)a3U(t,t')a3U{t',oo)] 

= tr[C/+(i,t>3C/(t,i')^3] , (20) 

where U{t, t') is an evolution operator defined in ([Ij. Ac- 
cording to the representation developed above, this func- 
tion is identical to the correlation function, K{t, t') — 
of the Majorana field, ^3(t), calculated with 
the use of action 53[^3, A] from Eq. (0). Importantly, in 
this formulation, the field ^3(t) is the representative of 
the Pauli matrix (T3. 

An observable of practical interest is the survival prob- 
ability, Pi.i{t,t') = \{i\U{t,t')\i)\^, which gives the 
probability of the spin to remain the same as it was in 
the beginning of evolution at time t. The survival prob- 
ability is related to the spin-spin correlation function via 
a simple identity: 



PUt^t') = ]^[l + K{t,t')]. 



Before proceeding, it is worthwhile to calculate the 
survival probability for A = and A2(<) = 0. In this 
limit the correlation function K(t,t') has a rather sim- 
ple form and can be calculated exactly. First we observe 
that one can perform Hubbard-Stratonovich transforma- 
tion in Eq. (jH]) by introducing an additional Majorana 



Pi(co,w) 
p,(o),-w) 



Pi,((0,w) 





2w 



-2 w 



p2 ((O.W) 



2(a),-w) 



FIG. 1: Diagrammatic representation for the matrix elements 
of G. Short [I3i{uj, ±w)] and long [I32{ljj, ±w)] arrows represent 
a; — 7- u; ± TO and lj — 5- cj ± 2w matrix elements respectively. 



field, 77, and then decouple the quadratic in Ai(i) term. 
In this new action. 



dt 



1 



766 



1 



r]fi + iA{t)^3 7] 



(22) 



the field Ai(i) plays the role of a "gauge field". Then 
we introduce variables = (^3 ± ri)/2 and eliminate 
the "gauge field," iAi{t) = W'^dtW with W{t) = 
exp[i J* dt' Ai{t')] = eyip[ix{t)] from the action by rescal- 
ing the fields, ^+ — M^^^. In this new fields action Eq. ([9]) 



acquires the form of a free field theory, iSy, 



Uodt 
transforms into 

K{t,t') 



25R 



but the correlation function now 



(|+(t)r(i'))e2^^t' -^''^it*') 

■t' 



sign[t - t'] cos 2i dt'Ai{t'] 



(23) 



This known expression can also be obtained directly from 
the original T-ordered exponent. 

If A ^ 0, it is convenient to perform calculations in fre- 
quency space. From Eq. ((TU]), it follows that the Fourier 
image of K{t,t') is nothing but the inverse of the op- 
erator, IC{uj,uj'). Therefore, the problem reduces to the 
calculation of this inverse matrix. 

Let us calculate K.{ui,ui') for the particular driving 
field in Eq. ([3]). Our strategy will be to firstly calcu- 
late G{uj,uj') from its definition (fTT|) . by inserting there 
the Fourier transformed driving field A-|-(a;) = A-{uj) = 
[e/2] S{uj) + A [S{uj -w) + S{uj + w)]. Then, integration 
over uJi yields 

G(a;,-w') = f3oiuJ,w)5{u ~Lu') (24) 
+ f5ki^,o'w)5{u! + akw — u)') 



fe=l,2;o-=± 



(21) where 



l3o{^,w) = e^nuj] + A^if[uj-w] + f[io + w]), 

(3i{u;,w) = eAif[u;] + f[u; + w]), 

I32{uj,w) = A^if[uj-w] + f[io + w]) (25) 

with f[x] = 2x/{x'^ — A^). Secondly, in order to find 
the inverse of IC~^ {i.e. K), we use the identity /C = 



5 



B(<o) Bo (CO) 

• = o + 



BoCm) 



Bo((a) B(co+w) B{Q)-w) Bo(Q)) 



B(tO) 

B(a) +2 w) B((B -2 w) 



B((B) 



B{(a) 
B„(a)) 



B{(b) 



FIG. 2: Dyson type equation for the diagonal elements of the 
matrix K, (full circles; see Fig. 3). Each of the five diagrams 
contributes to diagonal matrix element of K,. Thin arrows 
represent bare "hopping" elements /3i(aj,±w) and /32(<^,±™). 
Thick arrows represent corresponding dressed "hopping" ele- 
ments (see Fig. 4). 



^^p(l — /C^^)'^, where 1 stands for the identity matrix. 
This sum formally exist in the region |1 — /C^^| < 1. 
Outside of this region the inverse should be understood 
as the analytic continuation of the internal part. For 
our purpose it is convenient to develop a diagrammatic 
representation for the elements of the matrix }C~^ and 
its powers. We denote the powers of the diagonal 

elements of }C~^ as empty circle, elements /3i(a;, —w) 
and /?! (w, ui) as left and right arrows from uj to uj — w and 
w -\- w respectively. Similar arrows, but between uj and 
UJ ± 2w, will represent elements P2{uJ, ±w) (see Fig.l). 

Matrix elements of /C can be found by analyzing the 
structure of powers of 1 — . They have the following 
form 



/C(a;,a;') = B{u,w)5{lu - uj') (26) 

oo 

+ Tn{Lij,aw)5(ijj + anw — Lij'), 



where diagonal, B, and off-diagonal, r„, elements of K, 
will be determined below [see Eq. ([50]) ]. Importantly, 
Fourier transform of this expression yields the correlation 
function, K{t,t'). Diagonal and off-diagonal elements of 
JC satisfy certain functional relations. We identify these 
relations below and solve them exactly. 

Let us start with diagonal elements of K., namely, with 
the elements of the form B{uj, w)6{uj — lj'). These terms 
define the translational invariant part oi K{t, t'), which in 
time representation depends on the difference, t — t' , due 
to the presence of 6{uj — uj') in One can write Dyson 
type equation for the diagonal elements represented di- 
agrammatically in Fig. [2l where the full circles and the 
thick lines represent full series of diagrams presented in 
Figs. [3] and 2] respectively, and the empty circle, marked 
as Bq(uj,w), represents the full sum of bare diagonal el- 



(m-w) (to) (m + w) (m + 2 w ) ((o + 3w) 



B(co) 




FIG. 3: Diagrammatic representation for diagonal elements 
of the matrix /C. 



ements of the form [1 + 2uj ~ PoiuJ, w)] , leading to 



fe=0 



/?o(a;, w) - 2uj' 



(27) 



Functional relations pertinent to the diagrammatic series 
of Figs. [2] and Hread 

B{uj, w) = Bo{uj, w) - Bq{uj, w)B{uj, w) (28) 
X f3k{uJ,a-w)C{uj + aw, '-aw)B{uj + (Tkw,w). 

k=l,2;a=± 

Similarly, the off-diagonal terms in Eq. (|26l). that contain 
6{uj ^ uj' ±w) satisfy the relations following from Fig. 2] 

C{uj,w) = -Pi{uj,w) - l32iuj,w)C{uj + 2w,~w) (29) 
xB{uj + 2w, w) — C{uj, —w)j32{uj — w, w)B{uj — w, w). 

In Eqs. dlS]) and (HH), C{uj,±w) represents the fully 
dressed hopping matrix element of (1 — /C~^)~^, that cor- 
responds to the transition from a; to a; ± w. 

Finally, all remaining terms, Tn{uj, w), in Eq. (|26|) can 
be expressed via B[uj + nw, w) as follows 

n 

r„(w, w) = B{uj, w) Y[ C[uj ± (fc - l)w]B{uj ± kw, w). 



k=l 



(30) 



Solving Eq. (|28)) with respect to B{uj,w), one obtains 
it in terms of B{uj±w, w). This relation provides a possi- 
bility to generate continued fraction form of the solution 
for B{uj,w). Iterations in Eqs. P5|) and ([^ lead to the 
relations 

Co{uj,w) = -P^{uj,w), (31) 
C„+i(w,w) = -(ii{uj,w) - Cm{(^ + w,~w)l32{uj,w) X 

BjniuJ + 2w, w) - Crn{uj, -w)l32{uJ ~ W,w)Brn{uJ - W, w), 

B„i+i{uj, w)^^ = -2uj + l3o{uj, w) 

- ^ l3k{uj,crw)C„i{uj + aw,~crw)B„i{uj + aw^w). 

k=1.2-<y=± 
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C(m,w) P|((B,w) 



C (ra +2 w,-w) 



B((n+2w,w) 



,((B,w) 



C((B -w) 



B(a)-w,w) P2((a-w,w) 

FIG. 4: Dyson type equation for the non-diagonal (hopping 
elements of the matrix K.. 



While the solution of Eqs. (gH) and (HHI is simply de- 
fined by the infinite number of iterations B{uj,w) ~ 

which is the continued fraction representation of the di- 
agonal and off-diagonal matrix elements of /C(w, w')- 
alytical expressions for diagonal elements of K. obtained 
within two- and three-iteration approximation are pre- 
sented in Appendix B. 

Evaluating the matrix elements of K{uj) — 



J duj'K,{u 



which is the Fourier image of the correla- 



tion function K{t, 0), we find the spectrum of frequencies 
which contribute here. It is clear that all the matrix 
elements, r„(aj,±ui), of /C(aj,aj'), defined by Eq. 
will contribute to K{lij) substantially. It is also clear that 
due to the periodicity of B{u!,w) + J2'^=i^ni^,^w), 
as a function of w, which defines the Fourier image 
of the correlation function, K{t,t'), only a part of 
terms in the sum give essential contribution in the 
particular region of uj. Contribution of the tail becomes 
progressively smaller. We analytically calculate the 
first four elements of this sum, which gives the main 
contribution into the spin-spin correlation function in 
the region < u/A < 10 for the particular choices of the 
parameters A, w, a, and A. We have restricted ourselves 
within the fourth iteration level of the solution of 
Eqs. (l28l) for the same values of parameters. This means 
that we cut the exact continued fraction representation 
of IC{uj,uj') after four fractions, as we checked that five 
and more iterations do not affect the result for these 
parameters in the plotted range of uj/A. Comparison 
of our analytical expression for K{uj) with numerically 
evaluated solution of the corresponding Schrodinger 
equation are presented in Fig. [5] for various values of 
model parameters. Note that in these plots we took 
into account a finite relaxation rate, 7, in the Majorana 
fermion Green's function, which is needed to ensure 
the causality. More specifically we chose 7/A ^ 10~^ 
to ensure exponential decay of the Green's function at 
times t > 7^-'^. 



V. SUMMARY 

In summary, we have developed a new technique for 
studies of general non-equilibrium two-level systems. The 
technique is based on the mapping to a Majorana fermion 
field theory coupled to a scalar field. We have applied 



-1 ' 1 ' 1 ■ r 



3" 
i4 



< 

i4 




0.0 - 



FIG. 5: (Color online) Absolute value of the nonlinear di- 
electric response function plotted vs uj/A for various values 
of parameters. Full line is plotted from analytical expres- 
sions (see Eq. (|26p and below), while dotted line is obtained 
from exact numerics. Positions of the peaks determine the 
frequencies of the oscillatory Pi^i{t,t'). From top to bottom: 
a. w/A = 4, e/A=4, and A/A=2.05; b. w/A = 3, e/A = 4, 
and A/ A = 2.05; c. w/A = 4/3, e/A=l/3, and A/A = 1/6; 
d. w/A = 2/7, e/A=l/7, and A/A=l/28. Note that param- 
eters corresponding to figures c and d are out of the reach of 
the methods used in Ref. lid. 



the technique to study the dynamics of two-level sys- 
tems with driving fields given by Eq. ([3]). Our analyt- 
ical result for the nonlinear dielectric response function 
in energy space is shown to be in good agreement with 
the numerically evaluated solution of the time-dependent 
Schrodinger equation. We see that positions of the res- 
onance peaks in K{u}) are in agreement with the results 
of exact numerical simulations. 

Our technique allows generalization to dissipative two 
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level systems. We can consider decoherence being caused 
by the dissipative environment and also generated by dis- 
sipative elements in superconducting electronic circuits 
elements. One can extended our approach and include 
relaxation and dephasing times (Ti and T2) into consider- 
ation. However, these problems as well as the comparison 
of our approach to the well known comprehensive method 
based on the Bloch-Redfield equations^ is a subject of 
further studies. 
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VI. APPENDIX A 

Here we will present some of the basic properties of Ma- 
jorana fields, and will calculate the Green's function of a 
free non-interacting Majorana fermions. As fields, Majo- 
rana fermions are described by real Grassmann variables, 
^, with following propertie o^^'^^ 



(32) 



where curly brackets stand for an anticommutator. 

Integration rules over Grassmann variables are very 
simple. Namely 



(33) 



Using these rules of integration, one can prove that a 
Gaussian integral over Majorana fermions is expressed 
in terms of the determinant of the quadratic form. But 
contrary to ordinary fermionic integrals, it is equal to the 
square root of the determinant: 

/ 2?eexp{ J2 Cm^p-C.} = ^JlMA^] (34) 

This expression is correct both, for matrices and for dif- 
ferential operators. 



In order to calculate the Green's function, {^^i{t)£,i,{t')) , 
of non-interacting Majorana fermions, we introduce a 
generating functional 



V^a exp 



dt 



}. (35) 



From this expression it is clear that 



1 



(36) 



2o dri„{t)dr]f,{t' 

So we need to calculate Z{rj^). Calculation is straight- 
forward. One can find from Eq. (j35]) 



2(?7^ 



(37) 



X exp 



dt 



rl 



{^^, + 2i^^di^)dt[^^-2di\^) 



+ ^(ep-277,A"i)dt(e^ + 2drSM) 



where d^"^ = ^sign[t — t'] is the Green's function of 
the differential operator, dt = -^j- According to (p4t . 
the Gaussian integral over ± 2rj^dt^ yields Zo = 
^jDet[dtY, which is a C-number, and we obtain Z{ri^j) — 
Zq exp {— ?7^(i("^?7^} . Now, differentiating Z{r]fj^) twice 
with respect to 77^, and taking the limit Ty^ = 0, we re- 
produce the formula from the main text: 



(38) 



X exp 



Si_,^sign[ti - t2 



VII. APPENDIX B 

In this Appendix we present approximate expressions 
for solution of Eqs. (|28l) and ([2^ . which where obtained 
within two and three iterations {i.e. by cutting the con- 
tinued fraction representation after two and three frac- 
tions). 

Within two-iteration approximation we have for the 
diagonal element, Bi{ui,w), of matrix /C: 



Bi{uj,w)= (39) 

1 

n, ,!«„(',, ,„^ I , I3i(u,-w) h(t-^-w) ■ 

^ fjijyuj^ w j ^ _2(tj+«>)+/3o(w,w+u)) ^ -2{u}-w)+Mlo,lo-w) ~^ -2(u}+2w)+l}o{u} .,l^+2w) ' -2(tj-2u;)-|-/3o(tJ,cj-2ui) 

Similarly, for hopping element Ci{u)^w) we get from Eq. pil) 

r, ( N fl / ^ , /32(^,w)/3i {uj + 2w,-w) /?i(a;, -w)/32(^ - w, w) , 
C,{u:^w) = -P,iu;,w) + ^2iu + 2w)+Mc^ + 2w,w) + -2{u ~ w) + ^co - w)' ^^^^ 
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This yields the foUowing expression for diagonal elements of /C in a three-iteration approximation 
B2{uj,w) = 

1 

—2uj + I3q{uj,w) — Ci{uj,w)Bi{uj,uj + w) — Ci(a;, — — w) + I32{uj,w)Bi{lu,lu + 2w) + I32{l^, ~w)Bi(uj,uj — 2w) 

Finally, within the same level of approximation one has for hopping elements 6*2(0;, w): 

C2(a;, w) — — /3i(a;, w) — l32{uj, w)Ci{ll> + 2w, —w)Bo{uj + 2w, w) — Ci{ll>, —w)i32{uj — w, w)Bo{uj — w, w). (41) 

As we see from Eq. ([5Tt . one can continue this procedure to higher levels of iteration to analytically calculate the 
survival probability of the spin with an arbitrary accuracy. 
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